## Morrison 2009 IO

library(foreign) 
library(Amelia)
library(mi)
library(hot.deck)
library(mice)
library(mirt)
library(miceadds)

## Load original dataset
m2009 <- read.dta("M2009 IO Rep Data.dta")
head(m2009)
dim(m2009)

## Drop ID vars
m2009$countryname <- NULL

## Drop derived dummy vars
m2009$Dec70 <- m2009$Dec80 <- NULL

## How many variables? 62: no reduction necessary
dim(m2009)

## Imputation
## What is average percentage of missing data?
NAs <- function(x) {
    as.vector(apply(x, 2, function(x) length(which(is.na(x)))))
    }
NAs(m2009)
mean(NAs(m2009)/nrow(m2009))*100

## Thus: 37 imputations

## Note: already lags for nontaxtotalpc, grantspc, otherpc, soepc

set.seed(02138)
m2009.out <- amelia(m2009, m = 37, ts = "year", cs = "bankscode", idvars = "countrycode", polytime = 3, lags = c("politydv"), empri = 0.01*nrow(m2009))

write.amelia(obj=m2009.out, file.stem = "M2009 IO Imp Data", format = "dta", separate = FALSE)

## m = 5
set.seed(02138)
m2009.out.5 <- amelia(m2009, m = 5, ts = "year", cs = "bankscode", idvars = "countrycode", polytime = 3, lags = c("politydv"), empri = 0.01*nrow(m2009))

write.amelia(obj=m2009.out.5, file.stem = "M2009 IO Imp Data 5", format = "dta", separate = FALSE)

## m = % incomplete observations
complete <- complete.cases(m2009)
sum(!complete)/nrow(m2009)*100
## Thus: 94 imputations
set.seed(02138)
m2009.out.94 <- amelia(m2009, m = 94, ts = "year", cs = "bankscode", idvars = "countrycode", polytime = 3, lags = c("politydv"), empri = 0.01*nrow(m2009))

write.amelia(obj=m2009.out.94, file.stem = "M2009 IO Imp Data 94", format = "dta", separate = FALSE)

## MICE
m2009_mice <- mice(m2009, m = 37)
datlist <- mids2datlist(m2009_mice)
m2009_amelia <- list( "imputations"= datlist)
class(m2009_amelia) <- "amelia"
write.amelia(obj= m2009_amelia, file.stem = "M2009 IO Imp Data MICE", format = "dta", separate = FALSE)

## HD
m2009$countrycode <- NULL
head(m2009)
m2009.out.hd <- hot.deck(m2009, m = 37)
m2009.out.hd <- hd2amelia(m2009.out.hd)
write.amelia(obj= m2009.out.hd, file.stem = "M2009 IO Imp Data HD", format = "dta", separate = FALSE)
